MODULE Functions_libs

INTERFACE
  Subroutine TYPE1_DET(je,jht,jfix,jtype1) 
  Implicit NONE 
  Integer, INTENT(IN) :: je,jht,jfix
  Integer, INTENT(OUT) :: jtype1
  End Subroutine
END INTERFACE

 INTERFACE
  Subroutine TYPE2_DET(jdu,jdp,js,jtype2)
  Implicit NONE  
  Integer, INTENT(IN) :: jdu,jdp,js
  Integer, INTENT(OUT) :: jtype2
    End Subroutine
END INTERFACE

INTERFACE
  Subroutine TYPE3_DET(jh,jr,jtype3)
  Implicit NONE  
  Integer, INTENT(IN) :: jh,jr
  Integer, INTENT(OUT) :: jtype3
  End Subroutine
END INTERFACE

INTERFACE
  Subroutine TYPE4_DET(jmar,jos,jtype4)
  Implicit NONE  
  Integer, INTENT(IN) :: jmar,jos
  Integer, INTENT(OUT) :: jtype4
  End Subroutine
END INTERFACE


INTERFACE  
  Subroutine func_uti(jage,jtrpay,je,jh,jo,jmar,jtype2,CON,UTI)
 USE PARAMETERS
  Implicit NONE
Integer, INTENT(IN) ::  jage,je,jh,jo,jmar,jtrpay,jtype2
Real (Kind=dd), INTENT(IN) :: CON
Real (Kind=dd), INTENT(OUT) :: UTI
  End subroutine
END INTERFACE




INTERFACE  
  Subroutine func_uti_ret(je,jh,jmar,jtype2,CON,UTI)
 USE PARAMETERS
  Implicit NONE
Integer, INTENT(IN) ::  je, jh, jmar,jtype2
Real (Kind=dd), INTENT(IN) :: CON
Real (Kind=dd), INTENT(OUT) :: UTI
  End subroutine
END INTERFACE





INTERFACE  
  Subroutine Deterministic_z
  USE PARAMETERS
    Implicit NONE
  Integer:: x
  end subroutine Deterministic_z
END INTERFACE



INTERFACE
  function ran0(idum)
    USE PARAMETERS
    integer, INTENT(INOUT) :: idum
    Real (Kind=dd) :: ran0
  end function ran0   
END INTERFACE

INTERFACE
  function Tax_inc(income)
    USE PARAMETERS
Real (Kind=dd), INTENT(IN) :: income
Real (Kind=dd) :: Tax_inc
  end function Tax_inc  
END INTERFACE


INTERFACE
  function Bequest_val(beq_amount)
    USE PARAMETERS
Real (Kind=dd), INTENT(IN) :: beq_amount
Real (Kind=dd) :: Bequest_val
  end function Bequest_val
END INTERFACE


       
INTERFACE
  Subroutine maxxx(p,nr,nc,maxi)
    USE PARAMETERS
    integer, INTENT(IN) :: nr,nc
    Real (Kind=dd), INTENT(OUT) ::  maxi(nc)
    Real (Kind=dd), INTENT(IN) ::  p(nr,nc)
  end subroutine maxxx
END INTERFACE 


INTERFACE
  subroutine generate(idum,p,nr,rmaxi,indy)
    USE PARAMETERS
    Integer, INTENT(INOUT) :: idum
Integer, INTENT(IN) :: nr
    Integer, INTENT(OUT) :: indy
    Real  (Kind=dd), INTENT(IN) :: rmaxi,p(nr)
  end subroutine generate
END INTERFACE  

INTERFACE
Subroutine Moments(Array,n,ave,adev,sdev,var,skew,curt,cvar,logvar)
  USE PARAMETERS
Integer, Intent(In) :: n
Real(DD), Intent(In) :: Array(n)
Real(DD), Intent(Out) :: ave,adev,sdev,var,skew,curt,cvar,logvar
Integer:: j
Real(DD):: p,s,ep,temp
End Subroutine Moments
END INTERFACE 


END MODULE Functions_libs

!__________________________________________________________________________________________

 
  Subroutine TYPE2_DET(jdu,jdp,js,jtype2) 
  Implicit NONE 
  Integer, INTENT(IN) :: jdu,jdp,js
  Integer, INTENT(OUT) :: jtype2
if (jdu.eq.1) then 
jtype2=1
else 
jtype2=5
end if
if (jdp.eq.1) then 
jtype2=jtype2
else 
jtype2=jtype2+1
end if
if (js.eq.1) then  
jtype2=jtype2
else  
jtype2=jtype2+2
end if  
End subroutine


  Subroutine TYPE3_DET(jh,jr,jtype3)  
  Implicit NONE
  Integer, INTENT(IN) :: jh,jr
  Integer, INTENT(OUT) :: jtype3
  if (jh.EQ.1.AND.jr.EQ.1) then
jtype3=1
else if (jh.EQ.1.AND.jr.EQ.2) then
jtype3=2
else if (jh.EQ.2.AND.jr.EQ.1) then
jtype3=3
else if (jh.EQ.2.AND.jr.EQ.2) then
jtype3=4
else if (jh.EQ.3.AND.jr.EQ.1) then
jtype3=5
else if (jh.EQ.3.AND.jr.EQ.2) then
jtype3=6
end if
    End Subroutine
    
    
    
      Subroutine TYPE4_DET(jmar,jos,jtype4) 
  Implicit NONE 
  Integer, INTENT(IN) :: jmar, jos
  Integer, INTENT(OUT) :: jtype4
if (jmar.eq.1) then 
jtype4=1
else if ((jmar.eq.2).AND.(jos.eq.1)) then
jtype4=2
else if ((jmar.eq.2).AND.(jos.eq.2)) then
jtype4=3
end if  
End subroutine


    Subroutine TYPE1_DET(je,jht,jfix,jtype1)
  Implicit NONE  
  Integer, INTENT(IN) :: je,jht,jfix
  Integer, INTENT(OUT) :: jtype1
   if (je.eq.1.AND.jfix.eq.1.AND.jht.eq.1) then
     jtype1=1
     else if (je.eq.1.AND.jfix.eq.1.AND.jht.eq.2) then
     jtype1=2
     else if (je.eq.1.AND.jfix.eq.2.AND.jht.eq.1) then
     jtype1=3
     else if (je.eq.1.AND.jfix.eq.2.AND.jht.eq.2) then
     jtype1=4  
     else if (je.eq.1.AND.jfix.eq.3.AND.jht.eq.1) then
     jtype1=5
     else if (je.eq.1.AND.jfix.eq.3.AND.jht.eq.2) then
     jtype1=6
     
     else if (je.eq.2.AND.jfix.eq.1.AND.jht.eq.1) then
     jtype1=7
     else if (je.eq.2.AND.jfix.eq.1.AND.jht.eq.2) then
     jtype1=8
     else if (je.eq.2.AND.jfix.eq.2.AND.jht.eq.1) then
     jtype1=9
     else if (je.eq.2.AND.jfix.eq.2.AND.jht.eq.2) then
     jtype1=10
     else if (je.eq.2.AND.jfix.eq.3.AND.jht.eq.1) then
     jtype1=11
     else if (je.eq.2.AND.jfix.eq.3.AND.jht.eq.2) then
     jtype1=12
     
     
     else if (je.eq.3.AND.jfix.eq.1.AND.jht.eq.1) then
     jtype1=13
     else if (je.eq.3.AND.jfix.eq.1.AND.jht.eq.2) then
     jtype1=14
     else if (je.eq.3.AND.jfix.eq.2.AND.jht.eq.1) then
     jtype1=15
     else if (je.eq.3.AND.jfix.eq.2.AND.jht.eq.2) then
     jtype1=16
     else if (je.eq.3.AND.jfix.eq.3.AND.jht.eq.1) then
     jtype1=17
     else if (je.eq.3.AND.jfix.eq.3.AND.jht.eq.2) then
     jtype1=18
     
     end if
  
  End Subroutine

          


    
    !____________________________________________________________________________________________
Subroutine func_uti_ret(je,jh,jmar,jtype2,CON,UTI)
USE PARAMETERS
Implicit NONE
Integer, INTENT(IN) :: je,jh,jmar,jtype2
Real (Kind=dd), INTENT(IN) :: CON
Real (Kind=dd), INTENT(OUT) :: UTI

If (CON.le.0.0d0) then
  UTI=penalty
else if (jmar.eq.1) then
    UTI=   (      (CON**(alpha))   *  ((1.0d0 - phi_nw(jh,jtype2) )**(1.0d0-alpha))    )**(1.0d0-sigma)       /(1.0d0-sigma)
else if (jmar.eq.2) then
   UTI=   (((CON/equiv(je,n_ret))**(alpha)*( 1.0d0  - phi_nw(jh,jtype2)  )**(1.0d0-alpha) )**(1.0d0-sigma) ) /(1.0d0-sigma)

end if 

End subroutine

!____________________________________________________________________________________________
Subroutine func_uti(jage,jtrpay,je,jh,jo,jtype4,jtype2,CON,UTI)
USE PARAMETERS
Implicit NONE
Integer, INTENT(IN) :: jage,je,jh,jo,jtype4,jtrpay,jtype2
Real (Kind=dd), INTENT(IN) :: CON
Real (Kind=dd), INTENT(OUT) :: UTI

If (CON.le.0.0d0) then
  UTI=penalty
else if (jtype4.eq.1) then

    UTI=   (( (CON**(alpha))   *   ((leis(je,jo,jh,jtype2) )**(1.0d0-alpha))   )**(1.0d0-sigma) ) /(1.0d0-sigma)  -eta(jtrpay,jo)

else if (jtype4.eq.2) then

    UTI=   ((   ((CON/equiv(je,jage))**(alpha))  * ((leis(je,jo,jh,jtype2) + leis_mar(je,jo,1) )**(1.0d0-alpha)) )**(1.0d0-sigma) ) /(1.0d0-sigma)  -eta(jtrpay,jo)
    
    else if (jtype4.eq.3) then
        UTI=   ((   ((CON/equiv(je,jage))**(alpha))  * ((leis(je,jo,jh,jtype2) + leis_mar(je,jo,2) )**(1.0d0-alpha)) )**(1.0d0-sigma) ) /(1.0d0-sigma)  -eta(jtrpay,jo)

end if 

End subroutine





!_______________________________________________
!  function returns random number betweeen
!       0 and 1. Taken from Press et al (1994)
!_______________________________________________
function ran0(idum)
USE PARAMETERS
IMPLICIT NONE
integer, INTENT(INOUT) :: idum
Real (Kind=dd) :: ran0,am
INTEGER :: ia,im,iq,ir,mask 
parameter (ia=16807,im=2147483647,am=1./im,iq=127773,ir=2836,mask=123459876) 
integer k
idum=ieor(idum,mask)
k=idum/iq   
idum=ia*(idum-k*iq)-ir*k
if (idum .lt. 0) idum=idum+im 
	ran0=am*idum 
	idum=ieor(idum,mask) 
return 
end function ran0        

  function Tax_inc(jmar,income)
    USE PARAMETERS
Real (Kind=dd), INTENT(IN) :: income
Integer, INTENT(IN) :: jmar
Real (Kind=dd) :: Tax_inc
Tax_inc=max(0.0d0,  (( (income*5200.0d0) - tax_lambda(jmar) *( (income*5200.0d0)**(1.0d0-tax_tau(jmar))) )/5200 ))! for the purpose of calculating taxes, we convert to dollars, then back to model units.
end function Tax_inc 
    
 function Bequest_val(beq_amount)
    USE PARAMETERS
Real (Kind=dd), INTENT(IN) :: beq_amount
Real (Kind=dd) :: Bequest_val
Bequest_val= theta_beq * (( ( beq_amount + k_beq)**(1.0d0-sig_beq) ) / (1.0d0-sig_beq) ) 
  end function Bequest_val
          
Subroutine maxxx(p,nr,nc,maxi)
!       Given a shock today, finds the shock with maximum probability
!       among a matrix(nr,nc) with a probability distribution in each column
!       p is a matrix of dimension (nr,nc)
!       Output: maxi= maximum probability
USE PARAMETERS
IMPLICIT NONE
integer, INTENT(IN) :: nr,nc
Real (Kind=dd), INTENT(OUT) ::  maxi(nc)
Real (Kind=dd), INTENT(IN) ::  p(nr,nc)
integer jc,jr
Real (Kind=dd) ::  yy
!
Do 1 jc=1,nc
   yy=p(1,jc)
   Do 2 jr=2,nr
      if (yy.le.p(jr,jc)) then 
         yy=p(jr,jc)
      end if
 2    continue
   maxi(jc)=yy  
 1 continue
 return
 end subroutine maxxx
 
 
 
 !      Subroutine returns simulated current 
!      shock number(indy) for given probability distribution
!      DEFINITION OF VARIABLES
!      nr: number of shocks
!      p(nr): probability distribution
!      rmaxi: contains highest probability in the distribution (scalar)
!      Computed in subroutine maxxx
!      indy: number of simulated (drawn) shock


subroutine generate(idum,p,nr,rmaxi,indy)
USE PARAMETERS
USE Functions_Libs, ONLY : ran0
IMPLICIT NONE
Integer, INTENT(INOUT) :: idum
Integer, INTENT(IN) :: nr
Integer, INTENT(OUT) :: indy
Real  (Kind=dd), INTENT(IN) :: rmaxi,p(nr)
Integer ind
Real  (Kind=dd) :: ru,rx
ind=0 
1   if (ind.eq.0) then
       rx=ran0(idum)
       indy=int(nr*rx)+1
   	   if (indy.gt.nr) then 
	       print*, 'error indy', indy
	       print*, 'error rx', rx
	       print*, 'error nr', nr
	       print*, 'error idum', idum
           indy=nr
	   end if 
       ru=ran0(idum)
       if (ru.le.p(indy)/rmaxi) then
              ind=1
       end if
       goto 1
    end if
return
end subroutine generate


Subroutine Deterministic_z
USE PARAMETERS
Implicit NONE
Integer :: x

zd(:,:,:)= 0.0d0  ! educ, hc, health, past jo

!<=HS
!good health
do x=0,n_hc
   zd(1,x, 3)  = b0(1) + b1(1)*(dble(x)) + b2(1)*((dble(x))**2.0d0) + b3(1)*((dble(x))**3.0d0) + b5(1) 
end do

!average health
do x=0,n_hc
   zd(1,x, 2)  = b0(1) + b1(1)*(dble(x)) + b2(1)*((dble(x))**2.0d0) + b3(1)*((dble(x))**3.0d0) + b4(1)  
end do

!poor health
do x=0,n_hc
   zd(1,x, 1)  = b0(1) + b1(1)*(dble(x)) + b2(1)*((dble(x))**2.0d0) + b3(1)*((dble(x))**3.0d0)
end do


!Some college
!good health
do x=0,n_hc
   zd(2,x, 3)  = b0(2) + b1(2)*(dble(x)) + b2(2)*((dble(x))**2.0d0) + b3(2)*((dble(x))**3.0d0) + b5(2)
end do

!average health
do x=0,n_hc
   zd(2,x, 2)= b0(2) + b1(2)*(dble(x)) + b2(2)*((dble(x))**2.0d0) + b3(2)*((dble(x))**3.0d0) + b4(2)  
end do

!poor health
do x=0,n_hc
   zd(2,x, 1)= b0(2) + b1(2)*(dble(x)) + b2(2)*((dble(x))**2.0d0) + b3(2)*((dble(x))**3.0d0) 
end do

! college
!good health
do x=0,n_hc
   zd(3,x, 3)  = b0(3) + b1(3)*(dble(x)) + b2(3)*((dble(x))**2.0d0) + b3(3)*((dble(x))**3.0d0) + b5(3) 

end do

!average health
do x=0,n_hc
   zd(3,x, 2)=   b0(3) + b1(3)*(dble(x)) + b2(3)*((dble(x))**2.0d0) + b3(3)*((dble(x))**3.0d0) + b4(3) 

end do

!poor health
do x=0,n_hc
   zd(3,x, 1)=  b0(3) + b1(3)*(dble(x)) + b2(3)*((dble(x))**2.0d0) + b3(3)*((dble(x))**3.0d0) 

end do


end subroutine Deterministic_z



!**************************************************************************
Subroutine Moments(Array,n,ave,adev,sdev,var,skew,curt,cvar,logvar)
USE PARAMETERS
Integer, Intent(In) :: n
Real (Kind=dd), INTENT(IN) :: Array(n)
Real (Kind=dd), Intent(Out) :: ave,adev,sdev,var,skew,curt,cvar,logvar
	!Given an array Array(1:n) this subroutine returns the
	!mean ave, average deviation adev, standard deviation sdev,
	!variance var,skewness skew, kurtosis curt, coefficient of variation cvar,
	!and variance of the natural log logvar.
Integer:: j
Real (Kind=dd):: p,s,ep,temp


s=SUM(Array)
ave=s/n
adev=0.; var=0.; skew=0.; curt=0.; ep=0.
do j=1,n
	s=Array(j)-ave
	ep=ep+s
	adev=adev+abs(s)
	p=s*s
	var=var+p
	p=p*s
	skew=skew+p
	p=p*s
	curt=curt+p
Enddo
adev=adev/n
var=(var-ep**2/n)/(n-1)
sdev=sqrt(var)
cvar=sdev/ave
If (var .ne. 0) then
	skew=skew/(n*sdev**3)
	curt=curt/(n*var**2)-3.
Endif

ep=0.; logvar=0.
temp=log(ave)
do j=1,n
	s=log(Array(j))-temp
	ep=ep+s
	p=s*s
	logvar=logvar+p
Enddo
logvar=(logvar-ep**2/n)/(n-1)
End Subroutine Moments